#include "header.h"
#define PTOL         10.0
#define POLYL(x,y,z) (x+y+z)
#define POLY(x,y,z)  ( (x)*(x) + (y)*(y) + (z)*(z) + \
                       (x)*(y) + (y)*(z) + (z)*(x) + \
                       x + y + z )

/* ------------------------------------------------------------------------- */
void poly_test(struct node* p){
/* ------------------------------------------------------------------------- */
void poly_error(int, struct node*, double*, double*, double*, double*);
int i;
double ex, ey, ez, et;
double ex_max=0.0, ey_max=0.0, ez_max=0.0, et_max=0.0;
int ex_max_i=0, ey_max_i=0, ez_max_i=0, et_max_i=0;
FILE *fpt;

printf("Performing polynomial test, see POLY.OUT ...\n");
fpt = fopen("POLY.OUT", "w");
for(i=0; i<no_grid_points; i++){
    poly_error(i, p, &ex, &ey, &ez, &et);
    if(ex > ex_max){
        ex_max = ex;
        ex_max_i = i;
        }
    if(ey > ey_max){
        ey_max = ey;
        ey_max_i = i;
        }
    if(ez > ez_max){
        ez_max = ez;
        ez_max_i = i;
        }
    if(et > et_max){
        et_max = et;
        et_max_i = i;
        }
    fprintf(fpt, "%d %e %e %e %e\n", i, ex, ey, ez, et);
    if(ex>PTOL || ey>PTOL || ez>PTOL)
        printf("Node = %d: polynomial error HIGH, %f, %f, %f\n", i, ex, ey,
                ez);
    }
fclose(fpt);
printf("Maximum error in polynomial test:\n");
printf("x-derivative: %e at node %d\n", ex_max, ex_max_i);
printf("y-derivative: %e at node %d\n", ey_max, ey_max_i);
printf("z-derivative: %e at node %d\n", ez_max, ez_max_i);
printf("Total error : %e at node %d\n", et_max, et_max_i);
}

/* ------------------------------------------------------------------------- */
void poly_error(int pt, struct node* p, double *ex, double *ey, double *ez,
                double *et){
/* ------------------------------------------------------------------------- */
int i, neigh;
double x, y, z, poly, px=0.0, py=0.0, pz=0.0;
double dx=0.0, dy=0.0, dz=0.0, ds=0.0;

for(i=0; i<p[pt].total_conn; i++){
    neigh = p[pt].Conn[i];
    x     = p[neigh].x - p[pt].x;
    y     = p[neigh].y - p[pt].y;
    z     = p[neigh].z - p[pt].z;
    poly  = POLY(x, y, z);
    px += p[pt].X[i]*poly; 
    py += p[pt].Y[i]*poly; 
    pz += p[pt].Z[i]*poly; 
    dx = max(dx, fabs(x));
    dy = max(dy, fabs(y));
    dz = max(dz, fabs(z));
    ds = max(ds, sqrt(x*x+y*y+z*z));
    }

*ex = fabs(px - 1.0)/dx;
*ey = fabs(py - 1.0)/dy;
*ez = fabs(pz - 1.0)/dz;
*et = fabs((px-1.0)*(px-1.0) + (py-1.0)*(py-1.0) + (pz-1.0)*(pz-1.0));
*et = sqrt(*et)/ds;
}
